#include "stdafx.h"



/*template <>
static FILE& FILE::operator<<(FILE& f, Point& p)
{
        f << p.x << p.y;
        return f;
}

template <>
static FILE& FILE::operator>>(FILE& f, Point& p)
{
        if(f.GetVersion(1) < 2)...
}
*/
//***************************************************************



double Mesh::GetValue(Mesh &m, Point &p)
{
	int i;
	int f=1;
	double pp, s1, s2, s3, l1, l2, l3, s, u;

	Three_nodes *Tn;
	Four_nodes *Frn;
	Five_nodes *Fn;

	for (i=0; (i<kelem)&&f; i++)
	{
		if (m.NVTR[i].type == 4 || m.NVTR[i].type == 5)
		{
			if (p.x >= m.XY[m.NVTR[i].coord[0]].x && p.x <=  m.XY[m.NVTR[i].coord[1]].x && p.y >= m.XY[m.NVTR[i].coord[0]].y && p.y <=  m.XY[m.NVTR[i].coord[3]].y)
				f = 0;
		}
		else
		{
			l1 = sqrt((m.XY[m.NVTR[i].coord[0]].x-m.XY[m.NVTR[i].coord[1]].x)*(m.XY[m.NVTR[i].coord[0]].x-m.XY[m.NVTR[i].coord[1]].x)+(m.XY[m.NVTR[i].coord[0]].y-m.XY[m.NVTR[i].coord[1]].y)*(m.XY[m.NVTR[i].coord[0]].y-m.XY[m.NVTR[i].coord[1]].y));
			l2 = sqrt((m.XY[m.NVTR[i].coord[1]].x-m.XY[m.NVTR[i].coord[2]].x)*(m.XY[m.NVTR[i].coord[1]].x-m.XY[m.NVTR[i].coord[2]].x)+(m.XY[m.NVTR[i].coord[1]].y-m.XY[m.NVTR[i].coord[2]].y)*(m.XY[m.NVTR[i].coord[1]].y-m.XY[m.NVTR[i].coord[2]].y));
			l3 = sqrt((m.XY[m.NVTR[i].coord[0]].x-m.XY[m.NVTR[i].coord[2]].x)*(m.XY[m.NVTR[i].coord[0]].x-m.XY[m.NVTR[i].coord[2]].x)+(m.XY[m.NVTR[i].coord[0]].y-m.XY[m.NVTR[i].coord[2]].y)*(m.XY[m.NVTR[i].coord[0]].y-m.XY[m.NVTR[i].coord[2]].y));

			pp = (l1+l2+l3)/2.;

			s = sqrt(pp*(pp-l1)*(pp-l2)*(pp-l3));

			l1 = sqrt((m.XY[m.NVTR[i].coord[0]].x-m.XY[m.NVTR[i].coord[1]].x)*(m.XY[m.NVTR[i].coord[0]].x-m.XY[m.NVTR[i].coord[1]].x)+(m.XY[m.NVTR[i].coord[0]].y-m.XY[m.NVTR[i].coord[1]].y)*(m.XY[m.NVTR[i].coord[0]].y-m.XY[m.NVTR[i].coord[1]].y));
			l2 = sqrt((m.XY[m.NVTR[i].coord[1]].x-p.x)*(m.XY[m.NVTR[i].coord[1]].x-p.x)+(m.XY[m.NVTR[i].coord[1]].y-p.y)*(m.XY[m.NVTR[i].coord[1]].y-p.y));
			l3 = sqrt((m.XY[m.NVTR[i].coord[0]].x-p.x)*(m.XY[m.NVTR[i].coord[0]].x-p.x)+(m.XY[m.NVTR[i].coord[0]].y-p.y)*(m.XY[m.NVTR[i].coord[0]].y-p.y));
			
			pp = (l1+l2+l3)/2.;

			s1 = sqrt(pp*(pp-l1)*(pp-l2)*(pp-l3));

			l1 = sqrt((m.XY[m.NVTR[i].coord[2]].x-m.XY[m.NVTR[i].coord[1]].x)*(m.XY[m.NVTR[i].coord[2]].x-m.XY[m.NVTR[i].coord[1]].x)+(m.XY[m.NVTR[i].coord[2]].y-m.XY[m.NVTR[i].coord[1]].y)*(m.XY[m.NVTR[i].coord[2]].y-m.XY[m.NVTR[i].coord[1]].y));
			l2 = sqrt((m.XY[m.NVTR[i].coord[1]].x-p.x)*(m.XY[m.NVTR[i].coord[1]].x-p.x)+(m.XY[m.NVTR[i].coord[1]].y-p.y)*(m.XY[m.NVTR[i].coord[1]].y-p.y));
			l3 = sqrt((m.XY[m.NVTR[i].coord[2]].x-p.x)*(m.XY[m.NVTR[i].coord[2]].x-p.x)+(m.XY[m.NVTR[i].coord[2]].y-p.y)*(m.XY[m.NVTR[i].coord[2]].y-p.y));
			
			pp = (l1+l2+l3)/2.;

			s2 = sqrt(pp*(pp-l1)*(pp-l2)*(pp-l3));

			l1 = sqrt((m.XY[m.NVTR[i].coord[0]].x-m.XY[m.NVTR[i].coord[2]].x)*(m.XY[m.NVTR[i].coord[0]].x-m.XY[m.NVTR[i].coord[2]].x)+(m.XY[m.NVTR[i].coord[0]].y-m.XY[m.NVTR[i].coord[2]].y)*(m.XY[m.NVTR[i].coord[0]].y-m.XY[m.NVTR[i].coord[2]].y));
			l2 = sqrt((m.XY[m.NVTR[i].coord[2]].x-p.x)*(m.XY[m.NVTR[i].coord[2]].x-p.x)+(m.XY[m.NVTR[i].coord[2]].y-p.y)*(m.XY[m.NVTR[i].coord[2]].y-p.y));
			l3 = sqrt((m.XY[m.NVTR[i].coord[0]].x-p.x)*(m.XY[m.NVTR[i].coord[0]].x-p.x)+(m.XY[m.NVTR[i].coord[0]].y-p.y)*(m.XY[m.NVTR[i].coord[0]].y-p.y));
			
			pp = (l1+l2+l3)/2.;

			s3 = sqrt(pp*(pp-l1)*(pp-l2)*(pp-l3));

			if (fabs(s1+s2+s3-s)<1e-10) f=0;
		}

		u = m.NVTR[i].GetValueAnyPoint(m, p, i);
	}
}
/*double Three_nodes::GetValue(int Indx)
{
	return m->U[coord[Indx]];
}*/
double Three_nodes::GetValueAnyPoint(Mesh &m, Point *p, int Index)
{
	return U;
}

double Three_nodes::SetCoord(int &Array)
{
	coord[0] = Array[0];
	coord[1] = Array[1];
	coord[2] = Array[2];
}

/*double Four_nodes::GetValue(int Indx)
{
	return m.U[coord[Indx]];
}*/
double Four_nodes::GetValueAnyPoint(Mesh &m, Point *p, int Index)
{
	double res;
	
	double Psi1, Psi2, Psi3, Psi4;
	double Phi1, Phi2, Phi3, Phi4;
	
	Phi1 = m.XY[m.NVTR[Index].coord[1]]->Distance_x(p)/m.XY[m.NVTR[Index].coord[1]]->Distance_x(m.XY[m.NVTR[Index].coord[0]]);
	Phi2 = p->Distance_x(m.XY[m.NVTR[Index].coord[0]])/m.XY[m.NVTR[Index].coord[1]]->Distance_x(m.XY[m.NVTR[Index].coord[0]]);

	Phi3 = m.XY[m.NVTR[Index].coord[2]]->Distance_y(p)/m.XY[m.NVTR[Index].coord[2]]->Distance_y(m.XY[m.NVTR[Index].coord[0]]);
	Phi4 = p->Distance_y(m.XY[m.NVTR[Index].coord[0]])/m.XY[m.NVTR[Index].coord[2]]->Distance_y(m.XY[m.NVTR[Index].coord[0]]);
	

	Psi1 = Phi1*Phi3;
	Psi2 = Phi2*Phi3;
	
	Psi3 = Phi1*Phi4;
	Psi4 = Phi2*Phi4;

	res = m.U[m.NVTR[Index].coord[0]]*Psi1+m.U[m.NVTR[Index].coord[1]]*Psi2 + m.U[m.NVTR[Index].coord[2]]*Psi3 + m.U[m.NVTR[Index].coord[3]]*Psi4;

	return res;
}

double Four_nodes::SetCoord(int &Array)
{
	coord[0] = Array[0];
	coord[1] = Array[1];
	coord[2] = Array[2];
	coord[3] = Array[3];
}

/*double Five_nodes::GetValue(int Indx)
{
	return m.U[coord[Indx]]
}*/
double Five_nodes::GetValueAnyPoint(Mesh &m, Point *p, int Index)
{
		double res;
	
	double X1, X2, X3, X4, X5;
	double Y1, Y2, Y3, Y4, Y5;
	double Psi1, Psi2, Psi3, Psi4, Psi5;

	X1 =  m.XY[m.NVTR[Index].coord[1]]->Distance_x(p)/m.XY[m.NVTR[Index].coord[1]]->Distance_x(m.XY[m.NVTR[Index].coord[0]]);
	
	X2 =  p->Distance_x(m.XY[m.NVTR[Index].coord[0]])/m.XY[m.NVTR[Index].coord[1]]->Distance_x(m.XY[m.NVTR[Index].coord[0]]);
	
	if (p.x < m.XY[m.NVTR[Index].coord[4]].x)
			
		X3 = m.XY[m.NVTR[Index].coord[4]]->Distance_x(p)/m.XY[m.NVTR[Index].coord[4]]->Distance_x(m.XY[m.NVTR[Index].coord[0]]);	
	
	else X3 = 0;

	if (p.x > m.XY[m.NVTR[Index].coord[4]].x)
			
		X4 = p->Distance_x(m.XY[m.NVTR[Index].coord[4]])/m.XY[m.NVTR[Index].coord[1]]->Distance_x(m.XY[m.NVTR[Index].coord[4]]);	
	
	else X4 = 0;

	if (p.x < m.XY[m.NVTR[Index].coord[4]].x)
			
		X5 = p->Distance_x(m.XY[m.NVTR[Index].coord[0]])/m.XY[m.NVTR[Index].coord[4]]->Distance_x(m.XY[m.NVTR[Index].coord[0]]);	
	
	else 

		X5 = m.XY[m.NVTR[Index].coord[1]]->Distance_x(p)/m.XY[m.NVTR[Index].coord[1]]->Distance_x(m.XY[m.NVTR[Index].coord[4]]);

	
	Y2 =  m.XY[m.NVTR[Index].coord[1]]->Distance_y(p)/m.XY[m.NVTR[Index].coord[1]]->Distance_y(m.XY[m.NVTR[Index].coord[0]]);
	
	Y4 =  p->Distance_y(m.XY[m.NVTR[Index].coord[0]])/m.XY[m.NVTR[Index].coord[1]]->Distance_y(m.XY[m.NVTR[Index].coord[0]]);
	
	if (p.y < m.XY[m.NVTR[Index].coord[4]].y)
			
		Y1 = m.XY[m.NVTR[Index].coord[4]]->Distance_y(p)/m.XY[m.NVTR[Index].coord[4]]->Distance_y(m.XY[m.NVTR[Index].coord[0]]);	
	
	else Y1 = 0;

	if (p.y > m.XY[m.NVTR[Index].coord[4]].y)
			
		Y3 = p->Distance_y(m.XY[m.NVTR[Index].coord[4]])/m.XY[m.NVTR[Index].coord[1]]->Distance_y(m.XY[m.NVTR[Index].coord[4]]);	
	
	else Y3 = 0;

	if (p.y < m.XY[m.NVTR[Index].coord[4]].y)
			
		Y5 = p->Distance_y(m.XY[m.NVTR[Index].coord[0]])/m.XY[m.NVTR[Index].coord[4]]->Distance_y(m.XY[m.NVTR[Index].coord[0]]);	
	
	else 

		Y5 = m.XY[m.NVTR[Index].coord[1]]->Distance_y(p)/m.XY[m.NVTR[Index].coord[1]]->Distance_y(m.XY[m.NVTR[Index].coord[4]]);
	
	if ((m.XY[m.NVTR[Index].coord[4]].x!=m.XY[m.NVTR[Index].coord[0]].x)&&(m.XY[m.NVTR[Index].coord[4]].x!=m.XY[m.NVTR[Index].coord[1]].x))
		if (m.XY[m.NVTR[Index].coord[4]].y == m.XY[m.NVTR[Index].coord[0]].y)
		{
			Psi1 = X3*Y2;
			Psi2 = X4*Y2;
			Psi3 = X5*Y2;
			Psi4 = X1*Y4;
			Psi5 = X2*Y4;
		}
		else
		{
			Psi1 = X1*Y2;
			Psi2 = X2*Y2;
			Psi3 = X3*Y4;
			Psi4 = X2*Y3;
			Psi5 = X5*Y4;
		}
	else
		if (m.XY[m.NVTR[Index].coord[4]].x == m.XY[m.NVTR[Index].coord[0]].x)
		{
			Psi1 = X1*Y1;
			Psi2 = X2*Y2;
			Psi3 = X1*Y3;
			Psi4 = X2*Y4;
			Psi5 = X1*Y5;
		}
		else
		{
			Psi1 = X1*Y2;
			Psi2 = X2*Y1;
			Psi3 = X1*Y4;
			Psi4 = X2*Y3;
			Psi5 = X2*Y5;
		}

	res = m.U[m.NVTR[Index].coord[0]]*Psi1+m.U[m.NVTR[Index].coord[1]]*Psi2 + m.U[m.NVTR[Index].coord[2]]*Psi3 + m.U[m.NVTR[Index].coord[3]]*Psi4 + m.U[m.NVTR[Index].coord[4]]*Psi5;
		
	return res;
}

double Five_nodes::SetCoord(int &Array)
{
	coord[0] = Array[0];
	coord[1] = Array[1];
	coord[2] = Array[2];
	coord[3] = Array[3];
	coord[4] = Array[4];
}
/*
template <>
static FILE& FILE::operator<<(FILE& f, Element& e)
{
        int i;
		for (i=0; i<6; i++)	f << e.coord[i];
        return f;
}

template <>
static FILE& FILE::operator>>(FILE& f, Element& e)
{
        if(f.GetVersion(1) < 2)...
}
*/
//***************************************************************

/*template <>
static FILE& FILE::operator<<(FILE& f, Mesh& m)
{
        f << m.kuzlov << m.kelem;
        return f;
}

template <>
static FILE& FILE::operator>>(FILE& f, Point& m)
{
        if(f.GetVersion(1) < 2)...
}

template <>
static FILE& FILE::operator<<(FILE& f, double& d)
{
        f << d;
        return f;
}
*///*********************************************************************************************

main()
{

}